#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          //////////////////

  int Nf;
  double ainv;
  char dir[1024];

   int  set_p[4];
   char setpprint[16];
   char Oa2print[16];
   int iprop;
   int beta100;

   char momentum_order[1024];

          ////// Parsing command line


          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);

            if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);
            if(strncmp(argv[i],"-set_p",6)==0)  sscanf(argv[i],"-set_p=%d,%d,%d,%d",&set_p[3], &set_p[2], &set_p[1], &set_p[0]);
            if(strncmp(argv[i],"-iprop",6)==0)  sscanf(argv[i],"-iprop=%d",&iprop);
            
	    if(strncmp(argv[i],"-Oa2",4)==0)  sscanf(argv[i],"-Oa2=%s",Oa2print);


            if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-ainv",5)==0)  sscanf(argv[i],"-ainv=%lf",&ainv);

          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  whichProp(iprop,setpprint);
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum_Magro(momentum_order,set_p);
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  Equivalent_Momenta(SinP2);


	// IO

  char file[1024];
  char z_output[1024];
  char dir_qav[1024];
  sprintf(dir_qav,"%sa/%s/QAV/",dir,setpprint);
       //////////////////////////////////
        // Read Zq and LAMBDA
        /////////////////////////////

         double **Zq=malloc(sizeof(double*)*(Nclust+1));
	 for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double)*(Nmom_eq));

         double ***ZO=malloc(sizeof(double**)*(Nclust+1));
         for (int iclust=0; iclust<Nclust+1; iclust++) ZO[iclust]=malloc(sizeof(double*)*(Nmom_eq));
         for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) ZO[iclust][imom]=malloc(sizeof(double)*(5));

        double **aux_allmomenta=malloc(sizeof(double*)*(Nmom_eq));
       for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta[imom]=malloc(sizeof(double)*(Nclust+1));

        double *aux_clust=malloc(sizeof(double)*(Nclust+1));

    // Zq

    sprintf(file,"Zq_%s",Oa2print);
    for (int imom=0; imom<Nmom_eq; imom++)
    {
           read_input_onemomenta(dir_qav,file,aux_clust,imom);
           for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imom]=aux_clust[iclust];
    }


    // ZO 
   char gammaprint[8];

   for (int igamma=0; igamma<5; igamma++)
   {
         whichGammaProm(igamma,gammaprint);
	 sprintf(file,"Z%s_%s",gammaprint,Oa2print);
            for (int imom=0; imom<Nmom_eq; imom++)
            {
                read_input_onemomenta(dir_qav,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) ZO[iclust][imom][igamma]=aux_clust[iclust];
            }

    }
 

   ///////////////////////////////////////
   ////// COMPUTE ZRGI
   ///////////////////////////////////////


    double ***Zq_RIp_ainv=malloc(sizeof(double**)*(Nclust+1));
    for (int iclust=0; iclust<Nclust+1; iclust++) Zq_RIp_ainv[iclust]=malloc(sizeof(double*)*(Nmom_eq));
    for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) Zq_RIp_ainv[iclust][imom]=malloc(sizeof(double)*(1));

    double ***ZO_RIp_ainv=malloc(sizeof(double**)*(Nclust+1));
    for (int iclust=0; iclust<Nclust+1; iclust++) ZO_RIp_ainv[iclust]=malloc(sizeof(double*)*(Nmom_eq));
    for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom_eq; imom++) ZO_RIp_ainv[iclust][imom]=malloc(sizeof(double)*(5));


     
    for(int iclust=0; iclust<Nclust+1; iclust++)
    for (int imom=0; imom<Nmom_eq; imom++)
    {
		evolve_to_RIp_ainv(Zq_RIp_ainv[iclust][imom], ZO_RIp_ainv[iclust][imom], Zq[iclust][imom],ZO[iclust][imom],SinP2_eq[imom],Nf,ainv);
     }




        /////////////////
        // Write output
        ////////////////


      printf(" Writting output in  %s\n",dir_qav);

	double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom_eq));
	for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));
    ////////////////
    // Zq 
    ////////////////////
         sprintf(z_output,"%s/Z_ALLP/Zq_%s_RIp_ainv.dat",dir_qav,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zq_RIp_ainv[iclust][imom][0];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2_eq,Nmom_eq);

         sprintf(file,"Zq_%s_RIp_ainv",Oa2print);
         write_output_momentabymomenta_real(dir_qav,file,aux_allmomenta_real,SinP2_eq,Nmom_eq);
   ///////////////////////////////
   /// ZO R-PROM & GAMMA-PROM
   ///////////////////////////////////
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);

         sprintf(z_output,"%s/Z_ALLP/Z%s_%s_RIp_ainv.dat",dir_qav,gammaprint,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=ZO_RIp_ainv[iclust][imom][igamma];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2_eq,Nmom_eq);


          sprintf(file,"Z%s_%s_RIp_ainv",gammaprint,Oa2print);
         write_output_momentabymomenta_real(dir_qav,file,aux_allmomenta_real,SinP2_eq,Nmom_eq);

   }


	free(Zq); free(ZO);
	free(Zq_RIp_ainv); free(ZO_RIp_ainv);
	printf("ZRGI.c finished succesfully!");

    return 0;

}

